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Abstract. 

This paper extends an earlier quantum kinetics treatment for dilute, weakly- 
interacting, partially Bose-Einstein condensed gases, presented by the author elsewhere 
[J. Res. Natl. Inst. Stand. Technol. 101, 457 (1996)], by consistently treating 
the dynamics of the uncondensed atoms to the same level of approximation as the 
condensed atoms. Our method is based on a hierarchy of coupled equations of motion 
for the condensate mean field and fluctuations around this mean field, truncated 
to second order in the (effective two-body) interatomic potential, and with suitable 
decoupling approximations for higher order correlations. By applying perturbation 
theory in the Hartree-Fock-Bogoliubov basis, we re-derive the quantum kinetic theory 
of Walser et al. [Phys. Rev. A 59, 3878 (1999)], which further indicates the consistency 
of our treatment to the Kadanoff-Baym non-equilibrium Green's functions formalism 
for trapped gases. 
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1. Introduction 

The formulation of a consistent kinetic theory for the description of systems exhibiting 
Bose-Einstein condensation has been a very active area of research, since the very 
early days of the theoretical study of the behaviour of superfluid liquid helium [[[], |2|]. 
Such pioneering work for homogeneous systems has led to the establishment of the 
essential theoretical framework, which has since been extended to the case of trapped 
condensates. Due to the extreme weakness of their interactions, in comparison to those 
of liquid helium, the recently formed Bose-Einstein condensates in alkali gases are ideal 
systems for testing the validity of such theories. There are at the moment various non- 
equilibrium approaches to the dynamics of dilute, weakly-interacting, trapped Bose- 
Einstein condensates. Nonetheless, despite recent progress in the field, there exists 
to date no uniquely accepted non-equilibrium theory for the coupled dynamics of 
condensate and thermal cloud. On the one hand, we find the theories of Stoof [[J 
and Gardiner-Zoller et al. |4j], based respectively on a non-perturbative Fokker-Planck 
equation for the non-equilibrium dynamics of the gas, and a master equation for the 
many-body density matrix. These have been compared to experiments ||, yielding 
very good agreement between them and with existing data for the issue of onset of 
condensation || [7|], and they further appear to be in reasonable qualitative agreement 
with earlier work of Kagan, Svistunov and Shlyapnikov [0, |8|]. Other types of kinetic 
theories (except that of Castin and Dum assume the existence of a mean field 
for condensed atoms and deal perturbatively with more complex correlations around 
the mean field, based on suitable decoupling approximations, as first discussed by 
Kirkpatrick and Dorfmann [[!(]]. Such kinetic treatments have been recently discussed 
in the context of trapped gases by Proukakis and Burnett Hll], Zaremba-Nikuni-Griffin 



T2|| and Walser- Williams-Cooper-Holland |13], [14]]. In a recent paper, Wachter et al. 
have shown the equivalence of the non-equilibrium Green's function approach originally 
proposed by Kadanoff and Baym 0], as applied to inhomogeneous systems by Imamovic- 
Tomasovic and Griffin [[H], [17| , to the second-order gapless kinetic theory of Walser et 
al. [O, [l4j]- Central to this lies the fact that these theories yield the same second order 



damping rates, originally obtained in Beliaev's pioneering work |18| , and later extended 



by Popov [19], and more recently by Fedichev and Shlyapnikov [|20 



Based on a hierarchy of coupled equations of motion formalism for the condensate 



mean field and fluctuations around it, the author and Burnett |TTJ have established 
an equation for condensate dynamics which extends beyond Hartree-Fock-Bogoliubov 
theory by the inclusion of triplet correlations; these correlations are related to the source 
field generating the condensate in the Green's functions approach, as discussed, for 
example, in |T], [T2|, M. Discussion of these equations has led to the ab initio introduction 



of an effective interaction in them pifl , in a manner analogous to early diagrammatic 
work in the field [|18[ 19] . By focusing on static forms of these equations, Proukakis et al. 



22fl have further discussed gapless many-body theories, extending also beyond the bare t- 



matrix, establishing at the same time relations with other approaches [[23|, j24[. Extending 
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such work, Morgan |25j has formulated a gapless, number-conserving, self-consistent 
second order perturbation theory for the excitations of inhomogeneous condensates, 
which extends beyond the quasiparticle basis. Morgan has discussed at length the 
issues of gaplessness and infrared and ultraviolet divergences which have plagued the 
study of both homogeneous [18|, [HJ ^6], [27], |2B] and inhomogeneous condensates from 
their very beginning. 

In this paper, we extend the coupled equation of motion kinetic approach presented 
by the author elsewhere. In particular, in our preceeding work fTl] we discussed how 
one formulates equations beyond the simplest HFB theory, by discussing the equations 
of motion of correlations of up to three single-particle operators. Despite containing all 
the essential dynamics for the condensate mean field, that treatment only dealt with 
correlations of four single-particle operators in their mean-field approximation, hence 
not being able to generate the dynamics of the non-condensate in a consistent manner. 
To achieve this, one must additionally look into the dynamic off-equilibrium contribution 
of four-operator correlations, which drive such correlations away from their respective 
equilibrium mean values due to interactions with the surrounding particles. Doing this 
above the transition point, the author [2{| successfully obtained the classical Boltzmann 
equation, discussing at the same time how the actual interatomic potential becomes 
upgraded to an effective interaction. In this paper, we extend our earlier treatment by 
consistently treating, to second order in the effective (bare t-matrix) interaction, both 
the dynamics of the condensate and the non-condensate. This is based on the application 
of second order perturbation theory, starting from the usual Hartree-Fock-Bogoliubov 
(HFB) , or quasiparticle basis, in much the same manner as the theory of Morgan |f2H . 
Working in such a basis, we explicitly re-derive the second order theory of Walser et 
al. [|T^, [TJJ, which is of a different nature, formulated essentially as a perturbative 



expansion of a many-body density operator around a self-consistently adjusted Gaussian 
and its correlations. This equivalence was anticipated, since both theories rely on the 
assumption of a few slowly-varying quantities, given in our case, by the HFB order 



parameters. By further noting that the theory of Walser et al. [13|, |Tj] was recently 
shown |L5] to be equivalent to existing theories based on Green's functions |16|, [17J, we 
hence additionally establish a connection with such theories. Furthermore, since our 
perturbative treatment deals with the same hamiltonian as the equilibrium theory of 
Morgan ||25|| , we further believe that this paper establishes an indirect link between the 
static self-consistent approach of Morgan and the dynamic theory of Walser et al. We 
therefore consider this work as another step towards the formulation of a consistent 
theory for the study of both static and dynamic processes in finite temperature Bose- 
Einstein condensates. 

The layout of this paper is as follows. Firstly we present our formalism, giving 
explicit expressions for the treatment of the hamiltonian of the system and the 
corresponding energy functional. We then introduce, in Sec. Ill, the generalised matrix 
notation for the condensate and quasiparticle propagators and self energies, in the 
usual manner. In Sec. IV, we give explicitly the first and second order evolution of 
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these generalised propagators (in terms of the effective interatomic potential). Sec. V 
shows clearly the equivalence of our treatment to other kinetic theories (focusing on the 
recently proposed theory of Walser et al.), and we conclude this paper with some general 
remarks (Sec. VI). In the appendices, we explain how to deal with more complex normal 
and anomalous averages for both non-equilibrium contributions (via their respective 
equations of motion given in Appendix A) and their corresponding equilibrium values 
(for which we give the required decoupling approximations in Appendix B). 



2. Formalism 



We begin our treatment with the usual binary interaction Hamiltonian, 

H ^ ] ^rnd\d n 7) ^ ' V rsmn tt r (l s tt m CL n (1) 
rn rsmn 

which should incorporate most of the interesting physical processes occuring in dilute, 
weakly-interacting Bose-Einstein condensed gases. Here H = — (fr 2 V 2 )/(2m) + V trap (r) 
corresponds to the unperturbed hamiltonian in a harmonic trap, and V rsmn represents 
the symmetrised form of the interaction potential between a pair of particles, defined 
by Vrsmn = \ { (fs \V \mn) + (rs \V \nm) j . Here \i) = ipi(r) denotes a single-particle 
eigenstate of the trap, and the single-particle operators dj are related to the Bose field 
operator ^f(r,t) via $?(r,t) = J2i' l Pi( r )ci'i(t). Extending our earlier treatment regarding 
the ab initio introduction of an effective interaction [J21, 29\, Morgan |25) has shown that, 
as long as high-lying states are adiabatically eliminated, and one is only interested in 
states up to a certain cut-off [19|], the above hamiltonian can be equivalently written in 
terms of an effective re-summed two-body interaction. In fact, by numerically ensuring 
the independence of this effective interaction on the cut-off, this can be essentially 
replaced by the usual two-body t-matrix, often approximated by a local pseudopotential 
in three dimensions [|3(], 3JJ. All our subsequent expressions will be given in terms of 



V, bearing in mind that this essentially corresponds to such a re-summed two-body 
effective interaction. Our formalism relies on the existence of symmetry breaking, and 
hence we express the single-particle operators d<i £is 

h = (di) + (di - (di)) = Zi + ti (2) 

One can formulate an infinite hierarchy of coupled equations of motion for the 
condensate mean field and fluctuations around this field. To make this problem 
tractable, certain approximations are required. Firstly, we note that for the very dilute 
gases we are dealing with, most of the interesting physics should be already apparent by 
their second order expressions (in the effective re-summed interatomic interaction), since 
such expressions contain both energy shifts and irreversible damping processes. Our 
theory further requires consistent decoupling approximations (discussed in Appendix 
B). Due to the large number of atoms in existing condensates and the mean field 
potentials they generate, the single-particle eigenstates of the harmonic trap are not 
necessarily the best states for a perturbative expansion. Hence, starting from an HFB 
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basis, we incorporate all remaining effects as a perturbation. To this aim, we write our 
hamiltonian as 



H = H + H 1 + H 2 + H 3 + H 4 



where 



Hp — ^ rn z*z n + - v r . 



(3) 



(4) 



is the unpeturbed hamiltonian in a trap, additionally dressed by the condensate mean 
field, and the remaining operators are given by 

H\ ^ ] ^rn ^C^n ^ r C n -\- ^ ] Ksmn i^Z s Z m -\- (c^C m )^ Z ra -|- 2 r i^Z s Z m -\- (cjc m )^ C n 



"I - „ ^ ] Vrsmn (cJc^)c m Z ra -|- C^,Z S (c m C n ) 



(5) 



^2 ^ y 1 — ^™cjc n -|- ^ ] Vrsmn 2cJ ^Z g Z m -\- (cjc m )^ C n (cjc n )(cjc m ) 
rn rsmn 



„ ^ ] Vrsmn j^C^C^ {z m Z n -\- (c m C n )) -|- i^Z r Z g -\- (cjcj)^ C m C n {c^C^} (CmCnj (6) 



^3 ^ ] Ksmn ( 
rsmn 



(7) 



^ v cm n 



(8) 



The preceeding operator-dependent hamiltonians (H^) have been denoted by a prime, 
indicating that these are not the usual hamiltonians one would obtain from equation 
(H) upon writing a, = Zi + c\, but they have been modified by the application of the 
following mean field approximations (see e.g. |33], 34| ) 

cjcjc m — (cjc m )£g -|- (cjc m )cj -|- (cjcj)c m (9) 
cjc m c n — (cjc m )c n -|- (cjc n )c m -|- cj(c m c n ) (10) 

cJCgCj^jCn, — cj(cjx m )c n -|- cj(cjc n )c m -|- cj^cjxm)^ -|- C^C^CnjCm -\- (c^.C^ s ^)C m C n -\- C^,C^ s (^C m C n j 

- (2{cic n )(c\c m ) + (cicl)(c m c n )) (11) 

Such mean field approximations are conventionally performed in HFB treatments aimed 
at reducing the hamiltonian of the system to the usual quadratic form (^H + H[ + H'^J , 

ignoring completely contributions beyond the mean field labelled above as ( H 3 + H 4 



Ah 
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since, by definition, (\H' 3 + H'^j) = when taking averages over HFB eigenstates. Such 



a hamiltonian can be diagonalised by a transformation to quasiparticle operators [35 
which mix the shifted single-particle creation and annihilation operators c$ . In such a 
theory, the effects of thermal atoms are treated only in an approximate manner, and this 
is known to lead to an unphysical gap in the spectrum of elementary excitations in the 
homogeneous limit pi . This gap can be attributed to the inconsistent treatment 

of interactions between condensed and uncondensed components of the system and, more 
specifically, to the fact that the above approximation for three-particle operators is not 
physically justified ]^l], |25|, [36]]. For consistency we note here that the approximation 



on four operators arises from an extension of Wick's theorem |32[ (which is rigorous for 
equilibrium averages) to the non-equilibrium case. 

To avoid such complications, we will hence work with an "unperturbed" hamiltonian 
consisting of the usual quadratic HFB hamiltonian (^Ho + H[ + H'^J , but also keep the 

terms (h' 3 + H'^j such that our hamiltonian is still exact, as also done by Morgan ||25|| . 
As we will show explicitly in this paper, this implies that all first-order contributions 
(assuming no steady state anomalous correlations beyond the pair correlation) arise from 
the HFB hamiltonian (^H + H[ + H'^j , whereas the remaining parts of the hamiltonian 

(h' 3 + H'^J are solely responsible for all higher order collisional processes, damping 
rates, etc. In this paper, we will focus on the application of the above hamiltonian to 
the dynamical domain, somewhat similarly to the approach of Zaremba-Nikuni-Grifnn 
||12|| . We start our treatment from a suitably generalised version of our earlier kinetic 
equations, and gradually build upon them, which leads to the explicit re-derivation of 
the kinetic theory of Walser et al. [T3L WA . 

2.1. Some Definitions 

Applying perturbation theory on the HFB hamiltonian (^H + H[ + H'^J amounts to 
assuming that the correlations of up to two single-particle operators can be thought 
of as slowly-varying, whereas all other (higher order) correlations will be treated 
perturbatively (beyond their corresponding equilibrium contributions). We thus define 
the lowest order normal and anomalous averages of single-particle operators as 

Pa = (cfe> (12) 



K, 1 



(diCj) (13) 



These quantities, along with the condensate mean field Z{ form the set of order 
parameters of the HFB basis p2[ . This is analogous to the assumption made by Walser 



et al. [13] that such quantities form the set of relevant (master) variables for describing 
the evolution of the system on a coarse-grained timescale. For a binary hamiltonian, as 
already noted in [29|, a consistent treatment of the dynamic evolution of such quantities, 



to second order in the (effective) interatomic potential, requires an in-depth analysis of 
the evolution of all correlations, normal and anomalous, containing up to four single- 
particle shifted operators. As is well-known, higher order correlations containing an even 
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number of shifted single-particle operators do not vanish at equilibrium, but actually 
acquire a steady state value consistent with Wick's theorem. This states that, at thermal 
equilibrium, such higher order averages can be decomposed into products of lower order 
averages, which in our case correspond to the slowly-varying HFB order parameters 
z, p and k. To proceed further with our coupled equation of motion methodology, 
we must hence treat such quantities in two steps ||29|| : firstly we split them into their 
respective 'equilibrium' and 'non-equilibrium' components; subsequently, we decompose 
all equilibrium higher-order averages into products of HFB parameters, while at the 
same time perturbatively eliminating their rapidly evolving parts which drive them out 
of equilibrium. The non-equilibrium quantities we will be concerned with are defined 
by 

0~rsmn (cJ<3j(3 m C n ) ^cJcJc m C n ) (f^ r C^ s C m C n ^) (PmrPns PnrPms ^mn^rs") ( 1 '1 ) 

and 

£,rsmn (cJc s C m C n ) (c^C s C m C n } (& T C s C m C n ^) (Psr^mn Pmr^ns PnrPms) '"0 

This appearance of 'equilibrium' values is not encountered by triplet correlations (or 
higher order odd single-particle operator averages). These are assumed to vanish 
identically at equilibrium (due to our choice of the HFB basis as a good basis for 
describing the equilibrium properties of the system) and can therefore be treated 
perturbatively. In this case, we define p9| 



Armn (cJ<3 m C n ) (16) 
Tsrrara (CsC m C n ) (1^0 

2.2. Energy functional 

Having defined our notation, we can now explicitly write down the energy functional E 
of our system based on the hamiltonian (|I|). Its exact form is given by 

E=(H) = E$ FB + E^ (18) 



Here we have defined the HFB energy functional in the usual manner |32 

E HFB = ( (Ho + H[ + H' 2 ) ) = -rn (z*Z n + p nr ) 



~)~ „ ^ t Vrsmn \% r % s ^m %n ~\~ ^Pnr^ s Z m -\- < 2,p nr p m s ^rs^m^n Z r Z s Hmn ^ rs ^mh\Q) 
^ rsmn 

and the additional 'perturbative' beyond-HFB functional 

E^ = ((H 3 + H^j) — - ^ Vrsmn [^Knrs Z n + ^Kmn^* + 0" rsmn ] (20) 

^ rsmn 

Consideration of E^ FB leads to the well-known (reversible) first order expression for 
the dynamics of the HFB order parameters. The second term, E^ has no lowest order 
contribution in the HFB basis (since the HFB basis, by definition, does not allow for a 
non-vanishing steady state value for the triplets) and is hence conventionally ignored. 
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In this paper we demonstrate that it is precisely this second contribution to the 
energy functional which generates all second (and higher order) collisional contributions 
and should, thus, not be ignored P3| . In fact, we will show explicitly that the second 
order energy functional discussed by Walser et al. |L4| is nothing but our functional 
[e^ > fb + E^J when treated to second order, thus justifying its particular form. To 
obtain the desired higher order contributions in a perturbative expansion, we thus 
require the equations of motion of A and a, as given in Appendix A. 

Making use of these expressions and re-arranging, we obtain for the second order 
expression of the energy functional 

E {2) = H) E J2 V rsmnV pq lt 
rsmn pq lt 

[{.Pmp ^mp) {Pnq S n q) PtsPlr PmpPnq {Pts + 0~ts) {Plr + Sir)} 



X < 



+2 
+2 
+4 
+4 
+4 
+4 
+4 



{Pmp + S mp ) (ZqZn) PtsPlr ~ Pmp \ ZqZn) {Pts + $ts) {Plr + $l r ) 
\Pmp S m p) {pnq S n q) { Z s z t) Plr ~ PmpPnq { z * s Z t ) {Plr + Sir)] 
{pmp ~\~ S m p) ^nt^qsPlr Pmp^nt^q S {plr + Sir) 
{Pmp + S m p) KntKq S {z*Zl) Pmp^nt^qs { Z r Z l) 



^nt^qsPi 



Z p Z m 



KntK* {pir + Sir) 



{Pmp ~t~ S m p) { z n z t) ^qsPlr Pmp { z n z t) ^qs {Plr + Sir) 
{Pmp S m p) ^nt 



(21) 



Plr - PmpK-nt [ Z q Z s ) {Plr + S ir ) 

where 5^ is the usual Cronecker delta. 

We stress that the above expression only corresponds to the second order 
contribution of the energy functional E, and not to the exact energy functional given by 
equations (|T8|)-(f2"0|). In particular, we note that this second order expression additionally 
contains higher order averages such as 7, A, a and £ which play the crucial role of 
generating all higher order contributions to the energy functional (but do not modify 
the system evolution to second order, just as E^ has no effect on first order expressions 
for condensate / non-condensate evolution). 

In the above expression, we have defined the approximately energy- conserving 
matrix element 

%qit = J dt' e - i (^+^-^-^)(*-*') = 7c5{Au) - iP (-^j (22) 

where Aa; = {H>i + Q t — oj p — uj q ) and Ui denotes the energy of level i dressed by 
HFB mean field potentials (i.e. expectation value in the quasiparticle hamiltonian 
(^H Q + H[ + H^))- We note that the particular structure of the energy- conserving indices 
in the above expression for the energy functional (i.e. that they correspond precisely to 
the second matrix element V of E2) is a direct consequence of the Markov approximation 
which has been additionally imposed here (see section IV B). However, the Markov 
approximation only affects the indices appearing in the quasiparticle energies of the 
exponential in the integrand, and not the general structure of the energy functional. 

Apart from the first line of equation ( pl|) containing the usual Boltzmann 
scattering rates for normal (uncondensed) averages p, we find terms involving condensate 
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populations (z*z), as well as condensed (zz) and uncondensed (k) anomalous averages. 
The appealing feature of this energy functional (which, admittedly, is not apparent 
from its first order expression of equation ( pCf) ) is its explicit symmetry with respect to 
condensate and non-condensate contributions. For example, additionally to the usual 
'classical' Boltzmann rates (first line of equation (pip), we obtain all similar contributions 
in which one of the normal uncondensed averages p is replaced by (z*z). Even more 
appealing is the fact that this further holds when considering anomalous averages. 
Hence, in addition to terms ~ [(p + l)nn*p], E^ contains terms ~ [(p + l)n(zz)*p], 
[(p+ l)(zz)np\, [(p + 1)kk*(z*z)] and [{z*z)kk* p], etc. We note that it is precisely the 
choice of the HFB basis for our unperturbed hamiltonian which prohibits, in the final 
expressions, the appearance of multiple products of normal or anomalous condensate 
averages beyond their simplest forms (z*z), (zz) or (z*z*). This simplification would 
not arise if we were working with a simpler unperturbed hamiltonian, such as the 
one describing bare trap eigenstates, the Gross-Pitaevskii basis (dressed only by the 
condensate mean field) or the Hartree-Fock basis (additionally dressed by the normal 
non-condensate average). Finally, as already pointed out by Walser et al. |13[, although 



the in and out rates for the normal uncondensed component differ due to the process of 
bosonic enhancement (modifying p to (p+ 1)), the mean field (z*z) is never bosonically 
enhanced, and can thus be thought of as a classical field. 

3. Generalised matrix notation 

In order to present our formalism in the most general manner and establish a 
straighforward link with existing theories, we will henceforth work with generalised 
density matrices and hamiltonians, in a manner which clearly distinguishes between 
first order HFB and higher order perturbative results. We hence define the generalised 
condensate matrix R c (2n x 1) and the generalised quasiparticle density matrix R e 
(2n x 2n) by || 



*=U) *=U <A)J (23) 

where 1 is the unity matrix of the n-dimensional Fock space. These quantities can be 
thought of as the propagators for the condensed and uncondensed parts of the system, 
and are analogous to the corresponding generalised propagators conventionally employed 
in Green's functions formulations [[l], ||, [16|, |17], |4|, ^3], [37J. We note that, although R e 



includes all of the effects of a quasiparticle basis, these are explicitly given in terms of 
correlations of single-particle operators, with an explicit transformation to quasiparticle 
basis being beyond the scope of this paper (see e.g. |37|, |38|). 

Correspondingly, we define the generalised condensate and quasiparticle hamilto- 
nians as 

H - ( h(C) AW ^ H - ( h A \ (24) 

C ~l-(A( C ))* -{h&y J Ue -[- A * -h* 1 } 
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which can be thought of as the generalised self-energies for the above Green's functions. 

These generalised HFB matrix hamiltonians are already known from variational 
approaches, in which they can be readily obtained by variation of the energy functional 
Ejp FB with respect to the corresponding HFB propagators. In particular, one finds 



8E^ 



h v = raf = h % + T,m Vikij (z* k zi) = m\j) + 2 Em V Wj [z* k z t + p lk ] 

8E m (25) 

^3 = dK*.(t) = ^ij + ^ijkl (z k Zi) = J2kl Vijkl \ z k^l + «fcj] 



For consistency, and ease of subsequent comparison to the expressions of Walser 



et al. [14|, in this paper we have defined the generalised hamiltonians H c and H e in a 



slightly different manner from those of our earlier treatment JTI], This difference 

is only a matter of notation and arises from how the Pauli matrix a*- 3 -' = f ^, I 
3 \ -1 J 

is incorporated in our equations^]. We also point out the useful identities h*- = hji, 
Aij = Aji, pji = p*j and K k j = Kj k , whereas all higher order correlations are symmetric 
with respect to the interchange of indices labelling operators of the same 'type' (i.e. 
creation or annihilation c). 

4. Coupled evolution of condensate and quasiparticle propagators 

In this section we discuss the dynamic evolution of the condensate and quasiparticle 
propagators defined above, by employing the Heisenberg equation of motion (setting 
h = 1 for simplicity) 

i±(d) = ([6,H]) (26) 

for the mean value of a general operator O. 

4-1. Exact first order evolution 

To first order in the effective potential, the equations of motion for the condensed and 



uncondensed propagators acquire the following exact form [11, 



= H C R C + J (27) 
i^=(HR-RHA+K (28) 

Lib 

In each of the above equations, the first terms correspond to the HFB contributions, 
arising entirely from the hamiltonian (^Hq + H[ + H'^j . The static form of such equations 



can be readily derived variationally from E BFB |32| . Our approach extends beyond 



| We also note that, despite the different notation, our previous presentation contained a minor error, 
in that the left hand sides of equations (32) and (2.32) of [|ll| and |^9| respectively, should contain an 
additional factor of rj = a^). 
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such treatments, in that it further generates the 'non-equilibrium' matrices J and K for 
condensate and non-condensate evolution. These terms, arising solely from E^ 2 \ depend 
on averages of higher order than the HFB parameters, and are thus conventionally 
set to zero at equilibrium. In this work, these matrices are treated perturbatively, 
thus generating terms of higher order in the potential, which can still depend on the 
non-vanishing equilibrium quantities z, p and k. This shows that such contributions, 
and hence the beyond-HFB energy functional cannot be ignored even in static 



treatments, a point discussed clearly in the related work of Morgan p5| . To further 
establish a connection to kinetic treatments involving a damped nonlinear Schrodinger 
equation for the condensate mean field at finite temperatures, we note that it is precisely 
the presence of such additional contributions (labelled here by the kinetics matrix J) 
which lead to the damping term. In fact, by carefully identifying the physical origin 
of such a damping term (in a manner somewhat analogous to the present treatment), 
Zaremba et al. |L2[ and Bijlsma et al. || have respectively performed detailed studies 



of the hydrodynamic regime and the onset of condensation. It further appears that the 
kinetic matrix J is also taken into account in treatments based on linearised equations 
of motion, either explicitly |39] , or implicitly |34[ . 



Coming back to our formalism, the exact form of these kinetic matrices is given by 
L 

-L* 



J={ L „) (29) 



\-(N + N)* — (M — M*)* J v ; 



where 



smn 

Mji = ^ ] Vjsmn[o~ismn 4" ^imn z s ~\~ ^mis^n 4" ^nis^m] (32) 
smn 

Nji ^ Vjsmn\£,ismn ^sim^n ^sin^m r )imn^ , s \ (^'^) 
smn 

and iV represents the transpose of matrix N. Note that the above definitions follow 



the corresponding definitions of [fLl], £29], which they generalise by additional inclusion 
of quartic terms a and £ required for the consistent treatment of the non-condensate 
kinetics. 

4-2. Second order collisional integrals 

To proceed further and derive all relevant second order collisional integrals, we merely 
need to derive the corresponding equations of motion for the above kinetic matrices 
J and K. This requires the evolution of higher order normal and anomalous single- 
particle operator correlations, as given in Appendix A, and the application of consistent 
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decoupling approximations, given in Appendix B, for truncating the infinitely coupled 
equation of motion hierarchy. Our choice to apply perturbation theory starting from an 
HFB unperturbed basis, implies that the only effect of the hamiltonian (^Hq + H[ + H'^j , 
beyond yielding the first order evolution, is to define the renormalised unperturbed 
eigenenergies which appear in all higher-order expressions. In particular, we note 
that these eigenenergies dressed by HFB mean fields are implicit in the second-order 
collisional evolution of condensate / non-condensate propagators via the approximately 
energy-conserving matrix element V pq i t of equation (|22|) . To keep our subsequent 



notation compact, we henceforth suppress the free evolution of matrices J and K 
(see also Appendix A). Moreover, in giving their respective lowest order collisional 
evolution, we additionally impose the Markov approximation on the slowly-evolving 
HFB order parameters z, p and k appearing on the right hand sides of equations 
(|3~4^- (|35|) . This relies on the assumption that such quantities evolve significantly only 
over many complete collisional events, and thus decay much slower than all other 
interparticle correlations. We expect this to be a valid assumption for the dilute, weakly- 
interacting systems under consideration. (A detailed discussion of non-Markovian 
kinetic treatments which also deal explicitly with initial correlations in a system can 
be found in j|(J.) Within the above simplifications, we thus obtain 

.dj ( r 2 r z . 



dt 




(34) 



.dK 



dt 



[t pP + t k k* + i p ) 

T p k + T K (p + 1)* + I K 



h.c. 



(35) 



t p k + r K ( P + i)* + h 

[T pP + T k k* + I p ]* 

where h.c. denotes the hermitian conjugate and we have used a compact matrix notation, 
writing, for example, [(T p ) p\. { = J2i (Xp)jiPn^ an d similarly for the other contributions. 

We find that the rates T p and T K appearing in the above equation can be calculated 
variationally from the explicit second order expression for the energy functional of 
equation fl2"T|), as 



(T 



PJrl 



V dpi r J 



sq 



\ 9k* s J 

It is also interesting to note that 

( dE^ 



(36) 
(37) 



(I, 



Plrl 



\d(pir + Sl r ) J 

whereas we have not been able to obtain a correspondingly simple expression for I K 
For the condensate rates T <*) we find 



(3* 



( dE^ \ 



i 



hi z i 



>iq Z q 



dZi 

~dt 



(39) 
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where the subscript V 2 has been used to indicate second-order contributions to the 
evolution of the condensate mean field. Unlike the case of the condensate propagator, 



the presence of the non-vanishing matrices I p , I K in equation 



indicate that, in 



this case, one may not be able to define variationally a simple generalised hamiltonian 
including all rates (as was done earlier for the first order expressions). 

Below we define explicitly the form of all (n x n) matrices appearing in the above 
equations, in terms of their respective elements 



(T 
(T 



■ >il ^ Vj smn V pq lt 
smn pqlt 



[(Pmp 0~mp) {Pnq S n q) Pts PmpPnq {Pts + Sts)} 



+2 



{Pmp 0~rnp) ^nt^qs Pmp^nt^ 



qs 



)iq — 2 Vj smn Vpqlt {2 [(pmp ~\~ $mp) ^ntPls Pmp^nt {pis ~\~ Sis)]} 

smn pqlt 



(r P : 



2 ^ ' ^ t VjsmnVpqlt 
smn pqlt 



\{Pmp S m p) {pnq Snq ) Pts ~ PmpP I in (pts + ^ts)] 
")"2 {pmp 4" 0~mp) Pts Pmp i^ Z q Z n ) {pts + b~u 

4" [(Pmp "I - $mp) (Png &nq) { Z s Z t) PmpPnq { Z s Z t)\ 

+2 
+2 
+2 
+2 



(Pmp 4" 0~mp) ^nt^qs Pmp^nt^. 



{pmp 4" <5mp) { z n z t) Hq S Pmp { z n Z t) ^qs 
{Pmp ~\~ 0~mp) l^nt i^ Z q Z s ) Pmp^nt y Z q Z n 1 



(40) 
(41) 



(42) 



( r «)j 9 



2 ^ ] ^ ] VjsmnVpqlt 
smn pqlt 



2 [(Pmp 4" <5mp 



+2 
+2 
+2 



l^ntPls Pmp^nt {Pis + Sis)] 
{Pmp ~\~ 0~mp) l^nt { z *s Z l) Pmp^nt { z t Z l)\ 
( Z p Z rnj K ntPls ~ (^p^m) K nt {Pis + Sis) 
{Pmp ~\~ S~mp) { z n z t) Pis Pmp { z n z t ) {pis + Sis)] J 



PJji 



VjsmnVpqlt {Sli) 

smn pqlt 



PmpPnq {pts + S ts ) + 2 Z p Z m^j Pnq {pts + S ts ) + PmpPnq { Z s Z t) 
~\~^Pmp^nt^*qs ~l~ 2 {^ Z p Z nij ^nt^* qs ^Pmp { z n z t) ^* qs ^Pmp^nt ( Z q Z t 



{I*)ji 



2 VjsmnVpqlt (6u) 

smn pqlt 



(43) 



(44) 



x 



(45) 



2PmpK n t {Pis + Sis) + 2p mp K n t { z * s z l) + 2p mp {z n Z t ) {pi s + Si s ) 
+ 2 (zpZmj K n t {pig + Sis) + K-ml K ntKp S ^ K ml { z n z t) K p S + K-ml K nt ( Z p Z sJ 

Formal integration of the equations of motion (|34] ) -(|3~5| ) (taking into account their 
suppressed 'free' evolution in the HFB basis) and substitution into equations (p7D-(p8|) 
generates all second order collisional processes, as discussed in more detail in the next 
section. Before proceeding, we give at this point the most compact form of our equations 
to second order as 

(dRr\ ( dE 



dt 



dR r 



R r 



(46) 
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where the matrix / is defined as 
/ = 



I P h 



It I* 



h.c. 



14 
(47) 

(48) 



5. Link to other kinetic theories 

Having discussed the extent to which a variational methodology can be useful, let us 



now show explicitly how our equations reproduce those of Walser et al. |L3[], as written 
in generalised matrix form in their follow-up paper [|14|] . By thinking in terms of forward 
and backward scattering rates incorporating terms I p and I K , the equations of motion 
for the non-equilibrium matrix K can be straightforwardly re-expressed as 



.dK 
'~dt 



(rw p ) (r K + / K ) 



p « \ 

k* ( P * + i) J 

\f(p+l) K 

k K* p* 



+ h.c. (49) 



Following the notation of Walser et al. JHJ, we now define forward and backward 
collision operators as 



(50) 



and T > = — o\ (r < )*o"i. The corresponding generalised quasiparticle self-energies are 
= R e as defined in equation (23) and Rf = o~\ (Re)* &i- We thus obtain 



^ = - [(r<fl,< 



r > i?>) + h.c. 



(51) 



yielding for the quasiparticle propagator the following second order expression 

-iH e R> + (r<Rf - r>R>) + h.c. (52) 



dt 



Careful examination of the above equation shows clearly that our approach generates 
exactly the same forward and backward collision operators for the non-condensate 



propagator, as discussed in |14 



Although our methodology generates the forward and backward collision operators 
for the non-condensate propagator in a rather straightforward manner, the situation 
becomes somewhat more complicated when considering the evolution of the condensed 
component. The reason is that our methodology does not yield such expressions directly, 
but instead it provides us with the sum of forward and backward rates, via 



'dRc 
~dt 



' 8E 
dR r 



Rc. = -i 



dE 



(i) 

HFB 



+ 



dR c dR c 



R c = — iH c R c + YR C 



(53) 
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where 

" = ( -"[rli- ) (54) 

thus leaving their particular definition somewhat arbitrary. A natural definition 
corresponds to identifying these rates in terms of positive and negative contributions 
of r z (*) (denoted respectively by ±ve[r z (»)]), thus defining the forward and backward 
scattering rates for the condensate propagator as 

~ { -(+ve[V z *)y -(+ve[T z )y ) ^ 

and 

Y> = -a 1 (Y<)'a 1 =( + V %'\. , V"™) (56) 

V -(-ve[r z *\) -(-ve[T z \) J 

which clearly satisfy Y = (Y< — Y > ). 

In comparing our final expressions to those of Walser et al. [13, |14| , we note that 
our combined evolution equations agree entirely with the first formulation of the theory 
of Walser et al. fl3 |. In their subsequent work fll4| , where they have explicitly identified 
forward and backward rates as described in this section, their expressions for condensate 
rates differ from those we have given above in that their corresponding expression for 
Y z * of equation ( fH| ) is 

rr \ _oY^Y^T/ T> J 2 l(P™P + $rnp) K nt pi s ~ p mp K n t (pis + 0~ls)} 1 , - 

I 1 z*)iq - 1 2^ 2^ V jsmnVpqlt\ \ * „ 1 > [0 < ) 



smn pqlt 



^ml^nt^ps ^ml^nt^ps 



with the second line of the above expression, which identically cancels itself, being 
generated additionally to what appears in our expression of equation ([IT]). The 
appearance of such a contribution is physically appealing, since it implies that the 
collision processes occuring in condensate and non-condensate are of the same basic 
structure [JIJJ], indicating that processes T > can be generated from Y > by functional 
differentiation, as done in Green functions' techniques [jlj, 0. Although the net 
contribution of such a term is zero, its presence modifies the physical expressions of 
forward and backward collision integrals, i.e. it adds a particular collisional diagram to 
each. In our treatment, it is clear that such an extra 'self-cancelling' contribution cannot 
be generated in the first place. This is because the adiabatic elimination procedure 
yields for the second order expression for z terms proportional, at most, to correlations 
of five operators, whereas these additional contributions would require correlations of 
seven single-particle operators, and hence a treatment of three-particle interactions [20] . 
which extends beyond the binary hamiltonian employed here. It is therefore somewhat 
perplexing that such terms are found in the second formulation of the theory of Walser 
et al. |T3 |. However, if one is looking for a theory which satisfies such symmetries as 



mentioned above, due to specific physical considerations, we can indeed add and subtract 
these (or other similar) contributions to the forward and backward condensate collision 
integrals, obtaining exactly the same contributions as those discussed by Walser et al. 
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in |T4[ and by Wachter et al. JIR [38| . Such treatment will, however, not be a rigorous 
derivation, since its last stage will rely on additional physical arguments (e.g. gapless 
or conserving 16, |33|) arising from alternative approaches. 



We can now summarize the entire findings of our perturbative HFB treatment in 
the following systems of equations 



= H C R C + J 
= He R> - R>H\ + K 



i^ = -{Y<-Y>)R c 
m = - {(T<R< - T>R>) + h.c] 



(5f 



(59) 



Combining these results yields for the propagators of the system the following second 
order evolution 



4$* = -iH e R c + (Y<-Y>)R e 
t§- = -iH e R> + (T<R? - F>R>) + h.c. 



(60) 



6. Conclusions 



By applying second order perturbation theory from an HFB basis and making no 
further approximations on the binary interaction hamiltonian, we have generalised an 
earlier quantum kinetic approach to Bose-Einstein condensation presented by the author 
and Burnett elsewhere |TT| , by treating in a self-consistent manner the dynamics 
of both condensed and non-condensed components of the system. This involves a 
detailed consideration of the evolution of correlations of up to four shifted single- 
particle operators. Our methodology can be summarised as follows: Firstly we obtain 
the exact equations of motion for the condensate and non-condensate propagators in 
the binary interaction hamiltonian. These depend on matrices expressed in terms of 
higher order correlations, which can be approximately treated by separating off their 



'equilibrium' and 'non-equilibrium' contributions (somewhat analogously to [34 , |39l). 
The 'equilibrium' contributions are decomposed into lower order averages in terms of the 
HFB order parameters (thus yielding the usual first order evolution), whereas the non- 
equilibrium matrices contain all contributions of order higher than those obtained from 
the HFB energy functional. In this paper, we have focused on the lowest order expression 
for the evolution of such non-equilibrium matrices. Substitution of these results into 
the exact equations of motion for the condensate and non-condensate propagators yields 
their respective evolution to second order in the effective interatomic potentials and these 



are shown to be identical to those of Walser et al. [|K], [14| . Such an equivalence was 
of course anticipated, since they are both based on essentially the same underlying 
assumption of slowly- varying quantities. In the work of Walser et al. [[Uj, these 
quantities correspond to the 'master' variables determining the evolution of the gas 
on a coarse-grained timescale, whereas in the treatment presented here, these quantities 
are implicit in our choice of the HFB hamiltonian as the unperturbed hamiltonian for 
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applying perturbation theory. Our approach has provided the exact form of the energy 
functional (equations (|T8|)-(f20|)), which readily enables us to extend this treatement to 
higher orders. This may be useful, for example, to test a key underlying assumption 
of all such second order theories for dilute, weakly-interacting Bose-condensed gases, 
that the higher order terms (in the perturbative expansion) are significantly smaller. 
This should be valid sufficiently far from the critical region, when {kT /nU )(Vna 3 ) <C 1 



|T_9|, ^(J where a is the s-wave scattering length, n the condensate density and U Q the 
effective two-body interatomic potential (and can thus be justified by the fact that the 
problem has been formulated in terms of an unperturbed quasiparticle basis and that 
the interactions have been treated by an effective resummed bare t-matrix). 

Furthermore, our work is based on the same idea as the second order theory of 
excitations developed by Morgan |25] and recently applied numerically by the Oxford 
group [pi] , thus showing the link between these two approaches. We should also 
point out that our approach is formally connected to the detailed kinetic treatment 
of Zaremba-Nikuni-Griffin Jl2| (whose analysis focuses on the hydrodynamic regime). 
Moreover, in a recent paper the theory of Walser et al. |13|, [14]] has been shown 
to be equivalent to the Green's function methodology originally developed by Kadanoff 
and Baym and recently applied to the study of trapped condensates by Imamovic- 
Tomasovic and Griffin [TBI W7\. This latter approach is in turn formally related to the 
Schwinger-Keldysh method for treating non-equilibrium many-body systems, which has 
been employed by Stoof [|J. We stress that the work we have presented in this paper 
includes all relevant quasiparticle effects, despite being expressed in terms of single- 
particle operators; to generalise this, one could explicitly transform to quasiparticle 
operators. Carrying out such a transformation within the framework of the kinetic 
theory of Walser et al. [|T^, 13], Wachter has shown [03] that the collisional integrals 
simplify enormously, acquiring the usual form of Boltzmann-like factors, thus enabling 
a direct physical interpretation (which is rather concealed in the lengthy expressions 
in terms of single-particle operators). Similar work has been performed by Imamovic- 



Tomasovic and Griffin who recently described such a quasiparticle kinetic equation |37 
Here we should further point out that the static theory of Morgan p5 was formulated 
explicitly in terms of quasiparticles, and so where the related linear response discussions 
of Rusch et al. [39] and Giorgini pi]. The links of the current approach independently to 
those of Walser et al. (and hence Imamovic-Tomasovic and Griffin) and Morgan suggest 
that these time- dependent theories may correspond to some form of time-dependent 
generalisation of the gapless excitation theory of Morgan, although a more direct link 
remains to be established. Based on these considerations we hence believe that this 
paper represents yet another contribution towards the long-sought goal of a universally 
accepted kinetic theory for dilute, weakly-interacting systems exhibiting Bose-Einstein 
condensation. 
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Appendix A. Equations of motion for multiple shifted single-particle 
operator correlations 

In this appendix we give expressions for the equations of motion of correlations of 
up to four shifted single-particle operators. Our expressions are only explicit for the 
effect of the reduced hamiltonian (^H' 3 + H'^J , with the effect of the HFB hamiltonian 

(^H + H[ + H'2) being implicit in the renormalised eigenenergies cDj appearing below 
(and also in the approximately energy-conserving matrix elements V pq u of equation (22)). 
Explicit expressions in terms of the entire hamiltonian of the system H of equation ([!]) 
(including also the ones needed to derive the first order HFB-basis contributions to the 
coupled dynamics, i.e. first terms in each of equations fl2"T|)-(|2"8"D) can be found in jD], . 
We start by noting that the HFB quantities z, p, k will, to lowest order in the potential, 
have no contribution in the hamiltonian (h' 3 + H'^J , since the quantities found on their 
right-hand-sides acquire no equilibrium mean value within the HFB approximation. In 
particular, we find (making use of symmetries wite respect to summed indices) 

^"77 (^i) &i z i ^ t Vismn^smn (-^"-Q 
u smn 

i~7T (Pji) = {&j ~ &i) Pji + Vj smn {o~i smn + {2\* mig Z n + \imn z * s )] 

u smn 

5 j Vmnsi 0~mnsj i^^mjs z n ^jmn^s) (■^•2) 



i j, (^ji) (k'i ^j) ^ji 5 ] Vjsmn \£,simn {^^sim z n Ti'mn^s)] 

smn 

^ ] Vismn \£,sjmn (^^sjm z n 

+ 7,w£)l ( A -3) 

smn 

where the entire first order evolution of condensate / non-condensate propagators has 
now been summed up in the renormalised HFB eigenenergies uj{. 
For the non-HFB quantities, we find the following evolution 

i~77 (^smn) = {&m + &n ~ &a) ^smn + 52 ^mnlt (2pis z t) 
dt tt 

+ ^2 Vpnlt 2>PmpPls z t + ^^ml^ps^t + 2pis K mt Z p 
pit 
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^ dt ^ imn ^ 



+ V Vmlt 2p np pi s Z t + 2K n lK* s Z t + 2p U 
pit 

^ / Vpqls 2p m pP n qZi -\- 2 (yPmp^nl Pnp^ml) 
pql 



+ 2 X! Wmnlt { K il z t) + Vinlt ( K ml z t) + Vimlt (^nl z t)} 
It 

+ 2 Vpnlt Pip K ml z t + Pmp K il z t + ^il^mt Z p 
pit 

+ 2 Vpmlt PipKnl z t + PnpKil z t + ^il^nt Z p 
pit 

"i" 2 ^ ] ^piit Pmp^nl z t PnpK"ml z t 2fc m lfc n t Z p 
pit 



(A.4) 



(A.5) 



+ 51 Krmit (2piiPts) — (2pmpPnq) 



pq 



+ 2^2 V pnlt (PmpPliPts + Pls^mt K li + Pli K mt K p s ) 
pit 

+ 2^2 Vpmlt yPnpPliPts + Pls^nt^pi + Pli^nt^*p S 
pit 

2 ^ ] Vpqli {^PmqPnpPls 2p m pK n iK^ s 
pql 

2 ^ ] ^pg/s (PmqPnpPli 2p m pK n xK ( ^ 
pql 



(A.6) 



at 

+ X! Wmnlt (2pi s K it ) + Vi n i t (2pi s K mt ) + V im i t (2pi s K nt ) 
It 

+ 2 5^ V^j/t ^ml^it^ps 4" PmpPls^it + PipPls^mt 
pit 

+ 2 5^ V pm it ^nl^it^ps + PnpPls^it + PipPls^nt 
pit 

^ni^mt^ps "I - PnpPls^mt PmpPls^ntl 

pit 

2 ^ ] V^jgis [PrrapPrw?^ij "I - PmpPiq^nl PipPnq^ml] 
pql 



(A.7) 



We stress that equations (A.4)-( [^77| ) given above are not exact, even within our 
consistent decoupling approximations, because we have systematically ignored from their 
right hand sides contributions which depend on 'non-equilibrium' correlations of more 
than two shifted single-particle operators (i.e. A, 7, a, or £). Such contributions will 
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ultimately yield the evolution of the generalised condensate / quasiparticle propagators 
beyond second order in the potential, in the same manner that equations ( |A.1| )-( [AT3| ) 
will also yield contributions to third or higher orders. We believe this is fully equivalent 
to the statement made in Appendix A of Walser et al. [Of, that terms leading to 



higher order contributions are consistently neglected. We could, of course, have very 
straightforwardly included such terms in the equations of motion ( |A.4| )-( |A~7| ). However, 
since we do not deal with such terms any further, this would only introduce unnecessary 
complexity. 



Appendix B. Decoupling approximations 

As mentioned in the text, correlations of even numbers of shifted single-particle 
operators, can be thought of as consisting of two contributions [ 29fl , one part being 



due to the mean value at equilibrium when decomposed by Wick's theorem into 
products of lower averages (here referring to binary averages) and the remaining term 
yielding the rate of change of this quantity from its equilibrium value. For example 
{c\c\c m c n ) = ($ r c\c m c n )° + arsmn with the 'non-equilibrium' quantities (here a) treated 
perturbatively, by means of their respective equations of motion given in Appendix A. 

Since the HFB basis, by definition, does not allow for nonzero equilibrium values 
of triplet correlations, we note that, by analogy, any correlations of odd products of 
shifted single-particle operators will vanish to lowest order (since their final decomposed 
forms will always depend on such triplet products as discussed in [II, p9| ) and will yield 



non-vanishing contributions only at the next order in the potential, via their respective 
equations of motions discussed in Appendix A. The decoupling approximations we have 
employed for the 'equilibrium' quantities are as follows 

(c r C s C m C r ^j PmrPns PnrPms ^mn^rs (^"^) 

At equilibrium, this would arise directly from Wick's theorem, which is the main 
motivation for such a decoupling approximation here. The remaining decorrelations 
are performed analogously, namely 

(cldqClCt) = PqpKu + pi p K qt + p tp K q i (B.2) 
(CpCqClCt) = KpqKlt + KplK qt + K pt K q i (B.3) 



Pqp {PlrPts + PtrPls) + Pip {PqrPts + PtrPqs) + Ptp {PqrPls + PlrPqs) 
+ K q i [K^pts + K* s p tr + K* s p t p"j + K qt {K* pr pi s + K* s p tr + K* s pi p 

+K>lt \ K *prPqs + K psPqr + K rsPqp 



(B.4) 



'p Cq Cm CrClCt) 

Pmp (PrqKlt + Plq^rt + Ptq^rl) + Prp {Pmq^lt + Plq^mt + Ptq^ml) 
+Plp {P rq^mt Pmq^rt Ptq^mr) Ptp i^Prq^lm Pmq^rl Plq^mr) 
^rm^lt + Ki m K r t + KtmK r l) 



(B.5) 
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